MIKA: a multigrid-based program package for electronic structure calculations 
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A general real-space multigrid algorithm MIKA (Multigrid Instead of the K-spAce) for the self- 
consistent solution of the Kohn-Sham equations appearing in the state-of-the-art electronic-structure 
calculations is described. The most important part of the method is the multigrid solver for the 
Schrodinger equation. Our choice is the Rayleigh quotient multigrid method (RQMG), which applies 
directly to the minimization of the Rayleigh quotient on the finest level. Very coarse correction grids 
can be used, because there is in principle no need to represent the states on the coarse levels. The 
RQMG method is generalized for the simultaneous solution of all the states of the system using 
a penalty functional to keep the states orthogonal. Special care has been taken to optimize the 
iterations towards the self-consistency and to run the code in parallel computer architectures. The 
scheme has been implemented in multiple geometries. We show examples from electronic structure 
calculations employing nonlocal pseudopotentials and/or the jellium model. The RQMG solver is 
also applied for the calculation of positron states in solids. 



<N 



> 

in 
o 
m 
o 

<N 
O 

a 

O 

o 



X 



I. INTRODUCTION 

The goal of computational materials science and also 
that of modeling of nanoscale man-made structures is to 
calculate from first principles the various chemical and/or 
physical properties. This requires the solution of the elec- 
tronic (and ionic) structures of the system in question. 
The density-functional theory (DFT)tl makes a huge step 
towards this goal by casting the untractable problem of 
many interacting electrons to that of noninteracting par- 
ticles under the influence of an effective potential. How- 
ever, in order to apply DFT in practice one has to resort 
to approximations for electron exchange and correlation 
such as the local-density approximation (LDA) or the 
generalized-gradient approximation (GGA). Moreover, in 
the case of systems consisting of hundreds or more atoms 
it is still a challenge to solve numerically efficiently for 
the ensuing Kohn-Sham equations. 

We have developed a real-space multigrid method 
called MIKA (Multigrid Instead of the K-spAqe) for the 
numerical solution of the Kohn-Sham equationsB. In real- 
space methodsErO, the values of the wave-functions and 
potentials are presented using three-dimensional point 
grids, and the partial differential equations are di^ 
cretized using finite differences. Multigrid methodsD'H 
overcome the critical slowing-down (CSD) phenomenon 
occuring with basic real-space relaxation methods. Sev- 
eral approaches employing, the multigrid idea have ap- 
peared during recent yearsLTLj. 

From the different multigrid methods available for the 
solution of the Schrodinger equation, we have picked up 
the Rayleigh Quotient Multigrid (RQMG) method in- 
troduced by Mandel and McCormicldiJ r ^Xhis. r approach 
differs from full-approximation-storageE3l^'li3uj (FAS) 
methods, as well as from those methodsQ, where the 
eigenproblem is linearized. 

In the RQMG method the coarse grid relaxation passes 



are performed so that the Rayleigh quotient calculated 
on the fine grid will be minimized. In this way there is 
no requirement for the solution to be well represented 
on a coarse grid and the coarse grid representation prob- 
lem is avoided. Mandel and McCormicMlil introduced 
the method for the solution of the eigenpair correspond- 
ing to the lowest eigenvalue. We have generalized it to 
the simultaneous solution of a desired number of low- 
est eigenenergy states by developing a scheme which 
keeps the pjigenstates separated by the use of a penalty 
functional! 



II. NUMERICAL METHODS 



In our RQMG application the coarse grid relaxations 
are performed by the so-called coordinate relaxation 
method. It solves the discretized eigenproblem 



Hu = XBu 



by minimizing the Rayleigh quotient 

(u\H\u) 
(u\B\u)- 



(1) 



(2) 



Above, H and B are matrix operators chosen so that the 
Schrodinger equation discretized on a real-space point 
grid with spacing h is satisfied to a chosen order 0(h n ). 
In Eq. (H) u is a vector containing the wave function 
values at the grid points. In the relaxation method, the 
current estimate u is replaced by u' = u + ad, where 
the search vector d is simply chosen to be unity in one 
grid point and to vanish in all other points, and a is cho- 
sen to minimize the Rayleigh quotient. This leads to a 
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simple 1 quadratic equation for a. A complete coordi- 
nate relaxation pass is then obtained by performing the 
minimization at each point in turn and these passes can 
be repeated until the lowest state is found with desired 
accuracy. 



Non self-consistent solution 




Self-consistency iterations 




FIG. 1. Strategy of self-consistency iterations in MIKA. 
First, the wavefunctions are solved nonselfconsistently using 
the full multigrid method in the initial potential correspond- 
ing to the superposition of pseudoatoms. Then the effec- 
tive potential is updated (this is denoted by P in the fig- 
ure) . The potential update amounts to calculation of the new 
electron density, the solution of the Poisson equation and cal- 
culation of the new exchange correlation potential. Next the 
wave- functions are updated by one V- cycle. These two steps 
are repeated until self-consistency has been reached. 

Naturally, also the coordinate relaxation suffers from 
CSD because of the use of local information only in up- 
dating u in a certain point. In order to avoid it one 
applies the multigrid idea. In the multigrid scheme by 
Mandel and McCormickEil the crucial point is that coarse 
grid coordinate relaxation passes are performed so that 
the Rayleigh quotient calculated on the fine grid will be 



minimized. In this way there is no requirement for the 
solution to be well represented on a coarse grid. In prac- 
tice, a coarse grid search substitutes the fine grid solution 
by 

u'f = Uf + a//e c , (3) 

where the subscripts / and c stand for the fine and coarse 
grids, respectively, and // a prolongation operator in- 
terpolating the coarse grid vector to the fine grid. The 
Rayleigh quotient to be minimized is then 

{u f +od f c d c \H f \u f +aI f c d c ) _ 
{u f +al£ d c \Bf\uf +al£ d c ) 

(u f \H f u f )+2a(I c f H f u f \d c )+a 2 (d c \H c d c ) ( , 
(u f \B f u f )+2a(I c f B f u f \d c )+a 2 (d c \B c d c ) ' v 4 J 

The second form is obtained by relating the coarse grid 
operators, H c and B C1 with the fine grid ones, Hf and 
Bf, by the Galerkin condition 

H c = IjH f li; B c = I c f B f lf; If = (//) T . (5) 

The key point to note is that when HfUf and BfUf are 
provided from the fine grid to the coarse grid, the re- 
maining integrals can be calculated on the coarse grid 
itself. Thus one really applies coordinate relaxation on 
the coarse grids to minimize the fine level Rayleigh quo- 
tient. This is a major departure from the earlier meth- 
ods, which to some extent rely on the ability to represent 
the solution of some coarse grid equation on the coarse 
grid itself. Here, on the other hand, one can calculate 
the exact change in the Rayleigh quotient due to any 
coarse grid change, no matter how coarse the grid itself 
is. There is no equation whose solution would have to be 
represent able. 

In the MIKA package we have generalized the RQMG 
method to the simultaneous solution of several mutually 
orthogonal eigenpairs. The separation of the different 
states is divided into two or three subtasks. First, in or- 
der to make the coarse grid relaxations converge towards 
the desired state we apply a penalty functional scheme. 
Given the current approximations for the k lowest eigen- 
functions, the next lowest, (k + l)'th state is updated by 
minimizing the functional 

(u k + 1 \H\u k + 1 ) | |(^K + i)| 2 

(uk+i\B\uk+i) ~{ % (ui\ui) ■ (uk+i\uk+i)' 

The minimization of this functional is equivalent to im- 
posing the orthonormality constraints against the lower 
k states, when qi — > oo. By increasing the shifts qi any 
desired accuracy can be obtained, but in order to obtain 



*For the sake of simplicity, the wave- functions, and thus a, 
are here assumed real. We have implemented the complex 
case as well. 
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a computationally efficient algorithm a reasonable finite 
value should be used, for example 

Qi = (Afc+i - Ai) + Q, (7) 

where Q is a sufficiently large positive constant. In our 
test calculations Q is of the order of Q = 0.5 ... 2 Ha. 




FIG. 2. Electron density isosurface of the deep state local- 
ized at a neutral, ideal (no ion relaxation) vacancy in bulk 
Si. 

The substitution (||) is introduced in the functional (^) 
and the minimization with respect to a leads again to a 
quadratic equation. This time the coefficients contain 
terms due to the penalty part. 

While the penalty functional keeps the states sepa- 
rated on the coarse levels, we apply a simple relaxation 
method (Gauss-Seidel) on the finest level. The Gauss- 
Seidel method converges to the nearest eigenvalue, so 
ideally no additional orthogonalizations would be needed. 
In practice, however, we use Gramm- Schmidt orthogonal- 
izations and subspace rotationsH. However, the number 
of fine grid orthogonalizations remains quite plausible, 
for example, in comparison with the conjugate gradient 
search of eigenpairs employing only the finest gridtj. 

The Kohn-Sham equations have to be solved self- 
consistently, i.e. the wave functions solved from the 
single-particle equation determine via the density (so- 
lution of the Poisson equation and the calculation of 
the exchange-correlation potential) the effective poten- 
tial for which they should again be solved. To approach 
this self-consistency requires an optimized strategy so 
that numerical accuracy of the wave functions and the 
potential increase— in balance, enabling the most effi- 
cient convergencetj'B. Our strategy in MIKA for self- 
consistency iterations is illustrated in Fig. |l|. The Pois- 
son equation for the Coulomb potential is solved also by 
the multigrid method. 




FIG. 3. Left panel: Isosurface of the delocalized positron 
state in the perfect bulk SiCb- Right panel: Isosurface of the 
positron state localized at a Si vacancy in SiCb- 



III. EXAMPLES 

We have demonstratedi the performance of the MIKA 
scheme in calculating the electronic structures of small 
molecules and solid-state systems described by pseudopo- 
tentials. As a typical application Fig. || shows the elec- 
tron density of the so-called deep state localized at a 
neutral, ideal vacancy in bulk Si. It was shown, that the 
accuracy of 1 meV for the total energy was reached after 
three or four V-cycles, and that the amount of cpu-time 
needed was of the same order as when applying state- 
of-the-art plane- wave codes. We obtained an average 
convergence rate of approximately one decade per self- 
consistency iteration. This is,-etf the same order as those 
reported by Wang and Beekllj in their FAS scheme or 
by Kresse and FurthmullerEJj in their plane- wave scheme 
employing self-consistency iterations. The convergence 
rate of one decade per self-consistency iteration is bet- 
ter than that obtained by Ancilotto et aln in the FMG 
scheme and much better than the rate reached in the 
linearized multigrid scheme by Briggs et aln. 

We have applied the RQMG method also for the cal- 
culation of positron states in solids. Fig. || shows how 
the delocalized positron state in the perfect a-quartz is 
trapped in to a Si-vacancy. Positron states are a par- 
ticularly simple case for our method, because only the 
lowest energy wave function needs to be calculated in a 
given potential, so that no orthogonalizations or penalty 
functionals are needed. Moreover, in a simple scheme an 
electron density calculated without the influence of the 
positron can be used as the starting pointtll. However, 
even for the positron states the superior performance of 
the multigrid method in comparison with straightforward 
relaxation schemes is evident. For example, we have cal- 
culated the positron state in the Si-vacancy in bulk Si 
using a super cell containing 1727 atoms. The solution 
of the positron wave-function using the RQMG-method 
took less than a minute of cpu-time on a typical work 
station. To put this in proper context, J. E. Pask et afj 
report a similar calculation, based on the finite element 
method but without multigrid acceleration, for a super- 
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cell containing 4096 Cu atoms. The result converged 
within 1 ps took ' just 14.3 hr ' of CPU time. 




FIG. 4. Electronic structure of a four-carbon-atom chain 
between two jellium leads determined by the MIKA code. 
The figure shows an occupied sigma-type state localized at the 
atom chain. Periodic boundary conditions within the cylidri- 
cal symmetry have been employed along the (vertical) sym- 
metry axis and the (horizontal) radial directions. The centres 
of pseudo ions are denoted by grey circles and the positive 
jellium background by grey shadowing. Blue contours denote 
low electron densities whereas the density increases towards 
yellow and red. 

We have also applied the MIKA scheme in two- 
dimensional problems for quantum dots employing the 
current- spin- density functional theory (CSDFT), see 
Ref.tj. Moreover, we have implemented the RQMG- 
method in cylindrical coordinates enabling very efficient 
and accurate calculations for atomic chains, or systems 
which can be described using axisymmetric jellium mod- 
els. Fig. H shows a selected wavefunction of a system 
where a chain of four carbon atoms is sandwhiched be- 
tween two planar jellium leads. 

IV. SUMMARY AND OUTLOOK 

In the MIKA program package the-RQMG method in- 
troduced by Mandel and McCormicMlJ is generalized for 
the simultaneous solution of a desired number of low- 
est eigenenergy states. The approach can be viewed to 
belong to a third group of multigrid methods, in addi- 
tion to FAS and techniques where the eigenproblem is 
linearized. In principle, one can use arbitrarily coarse 
grids in RQMG, whereas in the other multigrid methods 
one has to be able to represent all the states also on the 
coarsest grid. 

We are convinced that our method will compete with 
the standard plane- wave methods for electronic structure 
calculations. However, some straightforward program- 
ming is still required. Implementation of the Hellmann- 
Feynman forces, required for the optimization of the ionic 
structures is under way. 



During the RQMG V-cycle, the states are all relaxed 
simultaneously and independently of each other. A paral- 
lelization over states would therefore be natural to imple- 
ment on a shared memory architecture. We have paral- 
lelized the MIKA codes over k-points, and over real-space 
domains. The domain decomposition is the appropriate 
method for distributed memory parallel computers. 

ACKNOWLEDGMENTS 

We acknowledge the contributions by Henri Saarikoski, 
Paula Havu, Esa Rasanen, Tero Hakala, and Sampsa Ri- 
ikonen in sharing their experience of the use of the MIKA 
package in different applications and preparing the fig- 
ures | (T.H.) and | (P.H.). T.T. acknowledges financial 
support by the Vilho, Yrjo and Kalle Vaisala foxunda- 
tion. This research has been supported by the Academy 
of Finland through its Centre of Excellence Programme 
(2000 - 2005). 



1 W. Kohn, Rev. Mod. Phys. 71, 1253 (1998). 

2 M. Heiskanen, T. Torsti, M. J. Puska, and R. M. Nieminen, 
Phys. Rev. B 63, 245106 (2001). 

3 T. L. Beck, Rev. Mod. Phys. 72, 1041 (2000). 

4 T. A. Arias, Rev. Mod. Phys. 71, 267 (1999). 

5 U. V. Waghmare, H. Kim, I. J. Park, N. Modine, P. 
Maragakis, and E. Kaxiras, Comp. Phys. Comm. 137, 341 
(2001). 

6 A. Brandt, Math. Comp. 31, 333 (1977). 

7 E. L. Briggs and D.J. Sullivan and J. Bernholc, Phys. Rev. 
B, 52, 5471 (1995); ibid., 54, 14362 (1996). 

8 F. Ancilotto, P. Blandin, and F. Toigo, Phys. Rev. B, 59, 
7868 (1999). 

9 J.-L. Fattebert, J. Comput. Phys. 149, 75 (1999). 

10 J. Wang and T. L. Beck J. Chem. Phys. 112, 9223 (2000). 

11 J. Mandel and S. McCormick, J. Comput. Phys. 80 (1989) 
442. 

12 A. Brandt, S. F. McCormick and J. W. Ruge, SIAM J. Sci. 
Stat. Comput. 4, 244 (1983). 

13 T. L. Beck and K. A. Iyer and M. P. Merrick, Int. J. Quan- 
tum Chem. 61, 341 (1997). 

14 S. Costiner and S. Ta'asan, Phys. Rev. E, 52, 1181 (1995). 

15 A. P. Seitsonen and M. J. Puska and R. M. Nieminen, Phys. 
Rev. B, 51, 14057 (1995). 

16 G. Kresse and J. Furthmiiller, Comput. Mat. Sci. 6, 15 
(1996). 

17 M. J. Puska and R. M. Nieminen, Rev. Mod. Phys. 66, 841 
(1994). 

18 J.E. Pask, B.M. Klein, PA. Sterne and C.Y. Fong, Comp. 
Phys. Comm. 135, 1, (2001). 

19 H. Saarikoski, M. J. Puska, and R. M. Nieminen, in this 
conference proceedings. 



4 



